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ABSTRACT 


Results are presented from a numerical investigation of the^ viscous 
transonic flow atout a stationary tody in free air. The geometry chosen 
was a symmetric NACA 64A010 airfoil at a freestream Mach number of 
0,8, a Reynolds number of 4 million based on chord, and angles of attack 
of 0- and 2 degrees. These conditions were such that, at 2 degrees 
incidence unsteady periodic motion was calculated along the aft portion 
of the airfoil and in its wake. This unsteady phenomenon has been 
experimentally observed for other airfoil geometries at transonic speeds. 
Although no unsteady measurements were made for the NACA 64A010 air- 
foil at these flow conditions, interpolated steady measurements of lift, 
drag, and surface static pressures compared favorably with corresponding 
computed time-averaged lift, drag, and surface static pressures. 
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SECTION' 1 


INTRODUCTION 

The principal ohjective of this research effort is to numerically 
calculate the viscous transonic flow field about a stationary lifting 
airfoil in free air. To our knowledge, this calculation represents 
the first-of-a-kind in transonic aerodynamics. 

Previous transonic lifting airfoil computations concerned steady 
inviscid flow with varying degrees of approximation to the Euler's 
eq_uations . ,^or example, Murman and Cole (Ref, l) , Krupp and Murman 
(Ref. 2), and more recently Murman, Bailey, and Johnson (Ref. 3) solved 
the transonic small disturbance sq.uation for flow past thin airfoils 
including imbedded shock waves . Transonic small disturbance theory 
solves the small perturbation equation (Ref. 4) with an approximate 
tangency condition at the airfoil surface, and under the assumption 
of isentropic flow, Jameson solved the transonic full potential 
equation (Ref. 5) using an exact tangency condition at the airfoil sur- 
face, and under the isentropic assumption (Ref, 6 ), Furthermore, 
Jameson's method permits computation of .imbedded shocks, 

Magnus and Yoshlhara (Ref. ?) solved the unsteady form of the 
Euler equations through an explicit finite difference scheme. A numer- 
ical solution of Euler's equations, including the possibility of embed- 
ded shocks (through shock capturing), embodies no other assumptions but 
the inviscid flow assumption. Coincident with the inviscid assumption, 
all three methods cited above employ a Kutta condition to determine the 
magnitude of the circulation for lifting bodies . 

For a lifting transonic airfoil with little or no separation in 
the neighborhood of the shock, combination of turbulent boundary 
layer theory (Refs. 8 and 9 ) with any of the above inviscid methods 
yields useful predictions of the shock location, lift coefficient, and 
drag coefficient, provided tliax a displacement thickness cjorrection 
is made to the airfoil shape. However, when the shock wave boundary 
layer interaction Induces considerable flow separation, 
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integrated inviscid-toundary layer theories are of dubious value 
in prediction of the shock location, lift, and drag of the system. 

Liepman and Roshko (Ref. lo) point out that transonic shock-wave 
boundary layer interactions involve an interplay between strength and 
position of shock waves on the body and boundary layer f character', i.e, , ' 

separation. The complete eq.uations of motion are thus needed to properly 
account for these nonlinear interactions . 

Numerical solution of the time-dependent Navier-Stokes equations, 
including turbulence and transition via models, yields predictions 
of the shock location,, lift coefficient, and drag coefficient of an air- 
foil under conditions of considerable flow separation. Thus, 
for thick transonic airfoils at moderate angles of attack, or for thin 
transonic airfoils at high angles of attack, numerical solution to the 
time -dependent Navier-Stokes equations is the only means for flow 
field prediction. 

A time -dependent Navier-Stokes calculation does not invoke a 
Kutta condition at the airfoil trailing edge, and Includes viscous 
effects ; hence , the possibility of vortex shedding , with attendent 
periodic flow is within the scope of this computation. For example, 

McDevitt, Levy, and Deiwert (Ref . ' ll) and Finke (Ref. 12) measured 
unsteady periodic flow about transonic airfoils for a limited range of 
Mach numbers and angles-of-attack. Finke ’s experiments on a symmetric 
NACA 63 ^ -012 airfoil at Mach ,70, Reynolds number 1.25 x 10^ based on 
chord, and two degrees incidence , show periodic flow at the 
airfoil trailing edge with a reduced frequency* K of 2.5 for the oscillating 
shock. The amplitide of shock oscillation is about one percent of the chore 
at two degrees angle of attack and increases to 10% of the chord at eight 
degrees incidence. 


■* The reduced frequency K is defined as the product of the angular 
frequency and the chord divided by the free stream velocity. 
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Finke (Ref. 12) also took interferogram pictures over one full 
period of shock oscillation for a q.uasi -elliptical NRl-profile, at a 
Mach numher of .?!, a Reynolds numher based on chord of 1 x 10 , and 
at 5 degrees angle of attack. These interferogram pictures, which 
represent lines of constant density, are shown in Figure 1. The 
pictures are depicted in intervals of 700/*s. In pictiire 3» the shock 
is positioned at 40 percent chord, in picture 5 the shock degenerates 
to a Mach wave and leading-edge separation determines the flow. Finally, 
picture 8 is similar to picture 2 indicating periodic flow. Shedding 
of vorticity is also indicated in the pictures of Figure 1, 

The geometry for this investigation was a symmetric NAGA 64A010 
airfoil at a freestream Mach numher of .80, a Reynolds numher of 4 x 
10^ based on chord, and at angles of attack of 0 and '2 degrees. At 
2 degrees incidence unsteady periodic motion was calculated along 
the aft portion of the airfoil and in its wake. Although no unsteady 
measurements were made for the NAGA 64A010 aiirfoil at these flow 
conditions, interpolated steady measurements of lift', drag, and 
surface static pressures compared favorably with corresponding 
computed time-averaged lift, drag, and surface static pressures. 

In order to solve for the unsteady viscous transonic flow field 
about a lifting body in free air, research was conducted in three 
principal areas. The are as follows; 

1. A computer code, called "STOKES", was developed which 
solved the Navier-Stokes equations with Multi-Regional Timesteps 
(MRT) and Vector-Do-Loops (VDL) . The MRT/VDL logic signifi- 
cantly reduced computational time. 

2. Turbulence and transition models were incorporated 
into the STOKES computer, code. 
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3. Bomdary conditions were generated along the perimeter 
of the iregion of calculation which simulated the lifting airfoil 
immersed in a frs6 airstreami 

A discussion of the STOKES computer code is presented in Section 
2, the turbulence and transition models embodied in the STOKES code are 
described in Section 3, boundary conditions are presented in Section 4, 
while results of the non-lifting (p^= O) and lifting “f) airfoil 

calculations are discussed in Sections 5 and 6, respectively. Section 
7 presents the conclusions reached in this research effort. 


SECTION 2 


THE STOKES CODE 

A computer code has "been developed, called "STOKES" , for calcu- 
lating viscous, compressihle, time-dependent flow fields about two- 
dimensional aerodynamic bodies. This code solves the Kavier-Stokes 
eauations including turbulence, transition, and free air boundary 
conditions along the perimeter of the domain of calculation. The 
finite difference eq.uations embodied in the STOKES code were originally 
developed by Trulio (Ref. 13) . These finite difference analogs of the 
eq,uations of motion are such that their self consistency property is 
maintained. That is, the finite difference equations for continuity, 
momenta, and internal energy imply an exact finite difference equation 
for total energy. The Trulio finite difference relations are second 
order accurate in space and first order accurate in time according to 
Taylor's series analysis. 

Turbulence and transition are discussed in Section 3> while the 
free air boundary conditions are described in Section 4. This Section 
is concerned with the numerical method employed to solve the finite 
difference analogs of the equations of motion embodied in the STOKES 
code. 

The principal innovation embodied in this computer code concerns 
computational time reduction. A method has been developed which produced 
computational time reduction factors between five and ten. For a non- 
lifting airfoil the flow field can be computed in less than 30 minutes 
on a GDC 76OO computer, while for a lifting airfoil’ the calculation 
requires less than two hours . . 

The STOKES computer code employs three principal time-reduction 
methods . 

1. Writing the program to proceed along quasi -streamlines 
- rather than the conventional processing along 
potential-like lines . 


2 . 


Incorporation of special assembly language subroutines 
to take maximal advantage of the "pipelining" capa- 
bility of the more advanced computers. 

3 . Setting up a procedure whereby the computational 

domain is automatically divided into several regions 
in which different timesteps govern numerical stability. 

2.1 Switched Computational Axis 

let us consider computation of the flow field about a two-dimen- 
sional body whose axis is parallel to the x coordinate axis , and 
immersed in an airstream moving along the x-axis . . The y-coordinate is 
considered the normal direction. The computation takes place on a finite 
difference mesh comprised of the intersection of q.uasi-streamlines , 
which are nearly parallel to the x-axis and extend upstream to downstream, 
and potential-like lines, which are nearly paralled to the y-axis, 
and initiate below the body and end above the body. 

Due to computer storage limitations, computations were previously 
conducted along the potential-like lines of the mesh, starting from 
the line farthest upstream of the body and ending at the line farthest 
downstream of the body. This mode of computation is referred to as 
unswitched. Usually the potential-like lines are shorter than the 
q.uasi -streamlines ; hence , they are comprised of fewer points . With 
points to a calculational line the do— loops in the program are 
executed fewer times, and there are correspondingly more interruptions 
for data transfer between central and peripheral memory banks. Further- 
more, most potential-like lines of the finite difference mesh pass through 
the body; hence, special branch point logic is required to account for 
the body aerodynamics... 
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The switched axis computer code perfoi^,- computations along 
the streamline -like lines of , the finite difference mesh rather than along 
the potential-like lines. For each calculational line the do-loops 
in the program are' executed a greater number Of times and' the number of 
data transfer interruptions are. reduced. Furthermore, the body 
geometry becomes two special streamlines of the mesh in the 
switched mode .of. computation.' Through most of . the' streamline 
computational sweeps there is no branching of logic. For the above 
reasons, the switched axis mode of computation is more efficient 
than the unswitched mode (Ref . 14 ) . 

2.2 The Vector Eo-Loop 

at NASA— Ames Research Center, Er. H. Eomax started devel- 
opment of GEC 7600 assembly language subroutines to replace operations 
done in regular FORTRAN do-loops. These subroutines were designated 
"Vector Eo-loops , " i.e., VEL*. The Vector-Eo-Loop takes advantage of 
the instruction stack of the GEC 76OO computer and employs the concept 
of pipeling in the computational seq.uence. In Reference l 4 , a simple 
vector do— loop is explained and an actual example is given showing the 
speed-up factor between a vector do-loop subroutine and the eq_uivalent 
FORTRAN program. 

The Vector-Bo-Loop substitutes computation during the waiting 
periods that normally, occur in FORTRAN do-loops. The case of the 
vector routine is a 'tight loop, with computations going on at different, 
stages of development, while the FORTRAN Bo-loop is essentially a series 
calculation. The' Vector-Eo-loop fits into the instruction stack 
and thus executes q.uickly. By making optimal use of the CHI registers, 
the Vector-Eo-Loop can do the work of a FORTRAN-Eo-Loop faster and 

more efficiently. . ... 

* The vector do-loops have also been referred to as the. SSPEEEY Codes. 
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2 . 3 Multi-Eegional Timesteps 

The Multi -Regional -Timestep (MET) logic divides the domain of 
computation into a set of concentric regions about the airfoil. In 
each concentric region a different timestep is employed, which is 
governed hy the stability of that particular region. The timestep 
generally increases by a factor of two as one moves from the inner 
concentric region towards the outer one. Thus, the inner regions are 
calculated frequently, while the outer regions are calculated in- 
frequently. 

The MET logic embodied in- the STOKES computer code is based on 
concepts developed by Magnus and Yoshihara (Refs. 7 and 15) and 
Trulio (Refs. l6 and 17). In solving for the inviscid field about an 
oscillating airfoil (Ref. 15) » Magnus and Yoshihara employed four con- 
centric mesh regions about the airfoil. Three are employed in the 
vicinity of the airfoil nose and the fourth is a coarse cartesian mesh 
which contains the other three and covers the remainder of the domain 
of calculation. The three in the nose region are a fine skewed mesh 
adjacent to the airfoil, a fine cartesian mesh containing the skewed 
mesh, and a medium cartesian mesh containing the other two. In their 
explicit, differencing scheme the allowable timestep is limited by 
stability criteria to an amount which is directly proportional to the 
spatial mesh increments. For each explicit timestep advancing the 
solution in the coarsest cartesian mesh of the system, 8 are made in 
the medium cartesian mesh, 64 are made in the fine cartesian mesh, and 
256 are made in the skewed elliptic mesh around the nose. Thus, 

Magnus and Yoshihara employed four, concentric regions in their compu- 
tation of flow about an oscillating airfoil. 

The multi-regional timestip formulation employed herein is patterned 
after that of Magnus and Yoshihara., A set of concentric regions is 
developed each of which extends from the upstream to the downstream 
boundary and contains the airfoil. In this foirmulation, which takes 


place in the switched axis computer code , the variahles of motion 
a.Te calculated for a given layer of zones as often as required by 
considerations of stability of that layer of zones and its immediate 
neighbors , rather than the former method of updating dependent variables 
for all zones • as often as required for the least stable zone of the 
finite difference mesh. Since the finite difference mesh is composed 
of streamline-like lines finely spaced in the neighborhood of a surface 
and coarsely spaced as one moves away from the surface, there are 
many different timestep levels in a typical aerodynamic problem. 

Therefore , a high percentage of the zones, i.e., those in the large 
timestep levels, will be calculated very infrequently and great savings 
is computational time should occur. 

The finite difference mesh is broken up into a series of regions 
comprised of a group of streamline -like lines. The regions are determined 
by stability for the least stable zone along a given streamline. The 
timestep levels monotonically decrease as one approaches the aero- 
dynamic body from below, and the timestep levels monotonically increase 
as one moves away from the top surface of the aerodynamic body. In the 
multi-regional timestep approach, one -dimensional variable timestep 
techniques are employed to sove a two-dimensional problem. Some of 
the techniques published by Trulio (Refs. i6,/7) are incorporated in 
the multi -regional timestep approach considered heirein. Particularly 
the concept of having adjacent concentric regions differ in timestep 
by a factor of two. Trulio has stated that under these conditions 
the conservation properties of the finite difference equations discussed 
in the introduction to Section 2 can be preserved after the solution 
is advanced in the coarsest region of the system, 

let us consider an airfoil finite difference mesh with four 
concentric regions, each differing by a factor of two in timestep, 
and let At be the minimum timestep, a "microstep". It is required 
to update the variables of motion on this mesh through the total time 


interval of 8 At, designated as a macro-step. The finest region will he 
updated through 8 increments in At, or through 8 " micros teps " , the 
next finest region will he updated through 4 timesteps of 2 At , the 
medium region will he updated through 2 timesteps of 4 At and the 
coarsest region through the macrostep. 

Complexity is minimized relative to the generalized variable 
timestep procedure of Trulio (Eefs. l4 and l6) , since all points along 
a given streamline belong to the same timestep level. Hence, since 
the streamline spacing controls the timesteps in most two-dimensional 
problems, the simple multi -regional timestep approach should yield 
most of the benefits of a more complex general two-dimensional variable 
timestep method. Development of the MET logic is presented xn much 
more detail in Eef erence 18 . 



SECTION 3 

TURBULENCE ANE TRANSITION MODELS 

An algebraic turbulence model, originally formulated by Baldwin 
and Rose (Ref . 19 )i is selected for the computation of turbulent flow 
over an airfoil at transonic speed. The model is, in effect, the 
mixing length theory to which relaxation along a streamline-like 
trajectory is incorporated. The simple algebraic model has been 
recently examined by Baldwin and Rose (Ref , 19 ) i (Ref. 20 ), 

and otheirs to be useful in shock boundary-layer interaction problems. 
All versions of available algebraic models are discussed. A criterion 
for boundary layer transition is also presented. 

3.1 Turbulence Modelling 

Numerical modelling of turbulence has become quite practical in 
the past decade with the advancement of high-speed computers. Though 
a universal model with wide range of applicability is far from reality , 
there is ample evidence that existing models have served well even in 
complex situations such as shock-wave boundary- layer interaction. All 
models of turbulence are supposed to be general in scope, and until 
recently, cross-comparisons between models (mainly studies done at 
NASA- Ames', Refs . 19 , 20 , and 21 ) are few . For transonic flow , there 
is no definitive conclusion as to the best turbulence model to employ. 

The usage of numerical modelis naturally bypasses the more funda- 
mental approach to turbulence studies via statistical theory, which 
might be at times academically pleasing but unrealistic in engineering 
applications. In general, turbulence modelling is divided into two 
categories: the algebraic models such as mixing length theory, and the 

transport models which are described by one^ or more differential equa- 
tions governing some quantity like turbulence energy, turbulence 
vorticity or shearing stress. The original work of Prandtl and its 
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sul3seq.uent ex^tension "by Geteci and Smith (CS) (Refs, 22 a.nd 23) > 

Rose (Ref. I9) , etc. are examples of the first class; the classical 
Kolmogorov model (Ref. 2 l^) (l9^2) and the Saffman model (1970) (Eef.25 ) 
fall into the latter category.' In adopting a transport model, one must 
solve, in addition to the "basic conservation laws, other differential 
equations from which tur"bulence stresses are determined. Transport .models try 
to depict the physics of turbulence tiansport, generation, dissipation 
and diffusion. In addition, some models (such as Saffman' s) show the 
correct analytical "behavior near the wall (as demanded "by the law of 
wall). The predictive capa"bilities for incompressi"ble "boundary layer 
flows "by those models are convincingly established. Turbulent flows 
in more than two spatial dimensions, including separation, compressi- 
bility, rotational effects, and containing boundary layers interacting 
with shock waves have not been subject to examination by those models*. 

In short, the transport models, as promising as they are, have yet to 
be * thoroughly tested 'by problems -more complex than plane boundary 
layer flows . 

In view of the existing complexities in the unsteady transonic 
flow problem, the desired economy in computation, and the added degree 
of complication in the nonlinear equations, we must see an alternative 
to the formulation by turbulence model equations. The alternative should 
be able to render a reasonably good description of the turbulent boundary 
layer development withou’fc a" di.sproportional amount of oomptft3’C.ion''il tii..e. 


* Wilcox (Ref. 26), applying Saffman 's model, has shown good results 
in the study of turbulent boundary separation and reattachment 
at moderate (2.96) Mach number. 


3.2 Algetiraic Turbulence Models 

The mixing length theory, originated by Prandtl (Ref. 27)1 pro- 
vides the foundation to all algebraic models. Modifications introduced 
by van Driest (Ref. 2S) , Cebeci and Smith (Refs. 22 and 23), and recently 
Baldwin and Rose (Ref. 19), Shang and Hankey (Ref. 29), and Deiwert 
(Ref, 20) all direct to improve the applicability of the model. Alge- 
braic models bypass the necessity of solving additional differential 
equations. From a computational standpoint, the eddy viscosity based 

on an algebraic model is post processed from mean-flow information. Our 
past application of the GS mixing length theory to internal flow 
problems in an impeller has shown good qualitative results (Ref. 30 ). 
Quantitative comparison is not possible due to the complete lack of 
experimental data. Hence, some version of an algebraic turbulence 
model is prefered to the more complex transport model. Despite the 
mixing-length common ingredient , there are variations in each indi- 
vidual formulation. The variations range from the unmodified theory 
to a relaxation model incorporating special treatment for the separated 
regions. The relaxation model was found significantly better than the 
unmodified algebraic model. According to Shang and Hankey (Ref. 29), 
it was significantly better than the Saffman's transport model for flow 
over a flat plate. Since separation on the transonic airfoil is a real 
possibility, incoirporation of the relaxation effect becomes quite 
desirable. For a detailed comparison of various formulations, we 
list them in the following table. 
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The formulation we shall adopt in the transonic flow problem is 
basically a hybrid relation primarily based on Rose’s relaxation 
model and the modification suggested by Delwert. Ingredients of 
the present algebraic model are blocked in heavy-lined rectangles 
in the preceeding table. The model for the turbulence stress 'XT^j 
can be svinunarized as follows. 
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( 1 ^ 


The eddy viscosity e is estimated by the mixing length theory 
which 'Subdivides the shear layer into an inner and an outer region. 
Inner region 


^ e = £ 2 




= k^yD 
k^ = 0.4 

y = normal distance from the nearest wall 
D = 1- exp (-y/A) 

A = 26v 


w 


/jM/e 


V = kinematic viscosity coefficient at the 


w 


nearest wall 


T^= shearing stress at the nearest wall 


Outer Region 


Selection 


^0 = 4 . 

= 0 . 095 ^ 

max 0 

= boundary layer thickness 


^eq. = 


C3) 


Relaxation along a streamline-like trajectory 


where 


\ = 


, if <fej < £ 


5i;, otherwise 
And 5 is a parameter defined along a streamline. 

Two major components, due to 'ieiwert ,, are introduced into Rome's 

|S)%f£5r in Pla« 
to avoid the complete vanishing of , e in a 


2Ui ^ dUr 



formulation. One is the adoption of 
of 

recirculating zone. Another one is the modification of the 

relaxation process in which e(^-A^) is used in place of a fixed e 

evaluated at' some reference station. Moreover, Deiwert found that 

relaxation over a streamline-like contour was more appropriate 

particularly for flows over a curved boundary, such as airfoil 

or turbine blade. Both modifications, indeed minor in nature, are 

convenient to implement with our computer code in which the 

scanning is done along streamline-like trajectories. The necessity 

Of incorporating the relaxation effect has been substantiated by 

Baldwin and Rose (Ref. ig) , Deiwert (Ref . 20) , Shang and Hankey (Ref . 29 ) . 

Its usefulness for flows In an unsteady transonic flow will be borne out 
our forthcoming computation. 
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3«3 Transition to Turbulence 


! 


Laminar flow at large Reynolds numbers becomes unstable, 
then the growth of disturbance in the boundary layer builds up 
until transition to turbulence occurs. The point of transition 
is strongly affected by the streamwise pressure gradient and the 
turbulence level of the free stream. To account for these factors, 
several empirical methods are available (for example, van Driest 
and Blumer, Crabtree, Granville, Smith and Gamberoni, van Ingen, 

Michel) . It is not possible to give a thorough comparison for 
those methods. In our blade-to—blado. computation for flows in an 
inpeller passage (Ref 50), toth Granville (-Ref 31),, and Michel's (Ref 32) formulations 
were examined. We found that Michel's simple algorithm provided 
a clear-cut prediction of transition point and it was extremely 
easy to implement. Since boundary layer transition is such a ' 


dubious subject in numerical computation, our guideline in the 
selection of a criterion is again "the simpler the better »" Unless 
future experiments contradict our selection, we shall adhere to 
Michel's criterion for the present application. The criterion 


gives a transition Reynolds number, (Re.), , based on the local 


Reynolds number Re . 


trcni. 


0 




The local Re^ can be estimated from the incompressible momentum 

(> - 

•>/ N ^K4X 


thickness 


e, = 


V. 




(SJ 
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transition to 


If the local Re^ is larger than (ReQ)trans, 
turbulent flow has taken place. Michel's criterion, resulted 
from correlation of experimental data, is supposingly valid for 
the range of Re between 0.1 x 10^ and 60 x 10 . Criterion of 
this nature signifies that transition to turbulence occurs at a 
point, rather than in a region, and relaminarization is not 
possible. A.M.O. Smith had compared Michel's algorithm against 
Granville's, Smith found the simple criterion of Michel quite 
satisfactory in the description of transition to turbulence. 



SECTION 4 


BOUNDARY CONDITIONS FOE SIMUIATION OF 
FREE AIR FLOW ABOUT AN AIRFOIL 


This research effort concerns flow ahout the symmetric NAGA 6UA010 
oil , which has a thickness to chord, ratio of IC^ « at two different 
angles of attack. At zero incidence the airfoil does not develop lift; 
hence, the circulation of the system is zero. For the 2 degrees incidence 
case the transonic field about this airfoil does contain circulation, 
since lift is being developed about the body. In this section boundary 
conditions are prescribed which simulate far field free air flow 
conditions for the non-lifting and lifting cases above. 


4.1 Boundaiy Conditions for Non-Lifting Airfoil 

Figure 2a shows a schematic of the domain of calculation of a 
symmetric airfoil having a chord G, a thickness T, and at zero degrees 
incidence. Let us consider inviscid, isentropic flow about this 
airfoil. .For the case of a thin airfoil, the flow is governed by the 
transonic small disturbance potential eauation (Ref. 


— 

where: 

XA04J (8) 

4y% 

<d =: ( 9 ) 
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( 12 ) 


) 

Let us define a scaled parameter y such that 
where ; 

<5 = 

ComMnation of Equations ( 7 ) and (l2) yields 


[k -i<+Q t ^ (13) 

where : 

( I - ( 14 ) 

4 ’s = (15) 

4)g = (16) 

€= 

Based on equations (lO), (ll) , (l5)f (l^)» 3-ad (l?), the streamwise 
and noirmal velocity components are, respectively: 

r: C4) y-H £ d)x ) (18) 

V- 


For the non-lifting symmetric airfoil of Figure 2a the effective 
far field potential solution to Equation (I 3 ) is a doublet. 



C 
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where D is the area of the airfoil section. Equation (20 ) indicates 
that the doublet potential strength diminishes inversely with the 
square of the transformed radius r = + Ky^)s from, the origin of 

the coordinate system, which is located at the 50 percent chord station 
(Figure 2a). Thus, the effects of the non-lifting airfoil on the 
perimeter of the system are small provided this perimeter is selected 
far from the airfoil. 

As a result of the above, the following boundary conditions were 
imposed on the perimeter of the domain of calculation shown in Figure 2a. 

(1) Along the upstream boundary the Murman-Cole inviscid 
solution (Ref. l) is prescribed for all time. 

(2) Frictionless flow is imposed along the lower horizontal 
lateral boundary of Figure 2a. 

(3) Frictionless flow is imposed on the dividing streamline 
upstream and downstream of the airfoil surface.! 

(4) No slip flow is imposed along the airfoil surface. 

( 5 ) A two-dimensional unsteady method of characteristics , 
including dissipation, is employed to determine the 
velocity field exiting from the downstream boundary (Ref. l4). 

At the upstream boundary the Murman-Gole solution is essentially 
the doublet of Equation (20). The lateral boundary is far enough _ 
away from the airfoil, so that any reasonable boundary condition will 
work including frictionless flow. F inally , at the downstream boundary 
the viscous equations of motion are solved, via the method of charact- 
eristics, to provide a very realistic model of the exiting flow. 

4.2 Boundary Conditions for lifting Airfoil 

A lifting airfoil develops circulation; hence, the boundary 
conditions enforced on the perimeter of the system greatly influence 
"the flow in the neighborhoo of the body. Let us consider the airfoil 
shown in Figure (2b) .at an angle of attack cL , and having the circu- 
lation r . The coordinate system (x, ) is located at the 50 percent 
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chord station, the, airfoil thickness is and its chord is G, 

Kirupp and Murman (Ref. 2) and Small (Ref. 33) have determined 
the far field potential solution to Equation (l2) for the lifting case. 
The resultant equation is as follows: 




0 - JiiH) Co£i± f... 




whersr; is the doublet potential given by Equation (20) and, 

©= 


( 22 ) 

(23) 


In Figure (2b) the point P is shown on the upper ; lateral boundary 
having the scaled radius r and the angle 0. The second term of 
Equation (2l) is directly dependent on the angle 0 and the circulation 
P ; hence, the effects of airfoil circulation are felt throughout the 
domain of calculation for the lifting case. Therefore, circulation 
must be properly accounted for far from the body to simulate free air 
boundary conditions. 

In order to preserve circulation . throughout the domain of calcu- 
lation, two sets of boundary conditions were investigated in this study. 
In the remainder of this section these boundary condition sets are 
described and their value ascertained. 

In the first set of boundary conditions, the velocity field on 
the upstream and lateral boundaries of the system shown in Figure (2b) 
was calculated at each timestep from the far field small disturbance 
Equations (18), (l9), (20), and (2l). The method of characteristics , 
which included dissipation, was utilized to compute velocities at the 
downstream boundary (Ref. 1^). 

To obtain a more accurate velocity field on the upstream and lateral 
boundaries of the system, the far field small disturbance equations 
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(18), (19)1 (20), (21) were solved in conjunction with the TSFOIl 
(Ref. 3) inviscid field, which served as the initial conditions 
for the lifting calculation. The TSFOIL velocity field was used 
to numerically evaluate the time-independent doublet contribution to 
the far field potential given by Eq^uation (2l), The numerically 
determined doublet contribution to the far field velocities , which 
resulted from the initial' conditions , was then saved and employed in 
Equation (2l) at later times. 

At each timestep of the calculation the lift force was computed 
about the airfoil and the local lift coefficient was determined, 

Based on the local lift coefficient the circulation about the airfoil, 
which was preserved along the perimeter of the system, was calculated 


from the following equation. 




r s Cl 


(24) 


Equation (24) is consistent with the small disturbance theory approx- 
imations. 

' The boundary conditions defined in this first set did not work. 
Vortices shed from the airfoil, whose rotation produced local velocities 
in the neighborhood of the upper lateral boundary (Figure (2b) ) 
opposite in sense to these prescribed at the lateral boundary from 
Equation (2l), Material could not flow out of the system; therefore, the 
pressure built up along the upper lateral boundary eventually wiping 
out the wake of the airfoil. 

As a result of the above, a second set of botindary conditions was 
prescribed in the far field. The velocity field on the upstream 
boundary was computed at; each timestep from the far field small dist- 
urbance Equations (l8), (19), (20), and (2l) in a manner identical 
to that described above. However, along the lateral boundaries of the 
system (.Figure 2b), far field small disturbance theory, via Equations 
(18), (19), (20), and (2l) , and the TSFOIL solution, were used 
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to compute only the external static pressure field along these 
boundaries . Based on the external far field theory pressures and the 
interior stresses along the lateral boundaries , the velocity field was 
computed from the equations of motion. Finally, the method of charact- 
0];7xstics, descirihed previously, was employed to compute velocities 
at the downstream boundary. 

The boundary conditions defined in the second set worked. During 
the shedding process the pressure boundary condition permitted material 
to flow out of the system at the upper lateral boundary. This per- 
mitted the shed vortices to pass through the system and resulted in 
a. reasonable description of., the flow field .for the wake of the airfoil. 
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SECTION 5 

RESULTS: 'OF NON-LIFTING CALCULATION 

Computations were first conducted for the NACA 64A010 airfoil 
at zero degrees angle of attack. Although the Reynolds number was 
4 X 10^, the computations were originally made with the turbulence 
model turned off; a computation is now in progress with the model 
turned on and with a finer mesh. The laminar problem was designated 

"Problem 101.0" 

The non-lifting airfoil numerical results for Problem 101.0 are 

presented in the following format; 

(1) Mesh used and Initial^ Conditions 

(2) Time-Histories of Pressure Coefficients 

(3) Computational Time Reduction Factor 

(4) Velocity Vectors of Steady Flow Field 

(5) Surface Pressure Distribution Comparisons 

(6) Mach Number Contours 

5.1 Mesh Used and Initial Conditions 

The mesh employed to solve Problem 101.0 is presented in Figure 3. 
This mesh is comprised of the intersection of 34 quasi-streamlines 
and 130 potential- like-lines. The streamline-like-line spacing is 
designed to provide about five points in the boundary layer downstream 

of the leeward airfoil shock. 

The timestep levels of the mesh of Figure 3 were investigated and 
the speed-up factor between a multi-regional timestep run and a con- 
stant timestep run was determined. Uniform freestream conditions were 
assumed in the numerical investigation. The speed-up factor is the ratio 
of the time it would take to do a macrostep (cover the total time 
increment of. the coarsest region) without the multi-regional time- 
steps to the time it would take to do the same computation with the 
multi-regional timesteps. Multi-regional timesteps and their speed-up 
factors are described in Reference l4. It was found that the minimum 


timestep for the mesh of Figure 3 was 1.33/e and the speed-up factor 
was 2 . 15 . Furthermore, it was foimd that the potential-like-line spacing 
in the vicinity of the leading edge of the airfoil controlled the time- 

step levels . 

A previous finite difference mesh was investigated, having the 
same number of mesh points, hut with finer potential-like-line spacing 
in the vicinity of the airfoil leading, edge. With finer potential-like- 
line spacing this finite difference mesh had a minimum timestep 
of .684^3 and a speed-up factor of I. 69 . Furthermore, this mesh is 
still severely limited hy the potential-like-line spacing in the 
vicinity of the airfoil leading edge, rather than the streamline-like- 

line spacing. 

Further study is warranted to increase the speed-up factor even 
further by spreading potential-like-lines near the airfoil leading edge. 

The initial field for the solution of Problem 101.0 was generated 
from the inviscid small disturbance solution from TSFOIL (Ref. 3 ). . 

The initial field generated from Murraan's method was interpolated 
onto the finite difference mesh of Figure 3 a computer code written 
especially for this purpose. 

A computer code, called "INTER", was developed to interpolate a 
flow field on a generalized finite difference mesh, termed the 
"unprimed points" onto another generalized finite difierence mesh, 
termed the "primed points". The sequence of events in the interpol- 
ation process are as follows; 

1. A primed point is located in the unprimed zone that contains 
it. 

2. A first order double Taylor's series expansion is employed 
to interpolate the data at the four corners of the unprimed 
quadrilateral onto the primed point. 

3. The velocity components, density, and specific internal 
energy 3-re interpolated. 
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The INTER computer code was employed to interpolate the 
inviscid flow field onto the 64A010 airfoil mesh of Figure 3. The 
inviscid field developed hy TSFOIL for the 64A010 airfoil was generated 
for zero angle of attack at Ifech .80. The inviscid field is developed 
on a rectangular mesh comprised of 50 horizontal lines and 89 vertical 
lines and considers the airfoil as a flat plate through use of the 
approximate tangency condition. The airfoil mesh (see Figure 3) is 
comprised of 3^ streamline-like-lines and I30 potential— like— lines . 

The interpolated pressure field on the 64A010 airfoil surface is 
presented in Figure 4. Measured pressure coefficients (Gross and 
Steinle , Ref 3^) at Reynolds number Al- x 10^, are also included for 
comparison purposes. It is seen from Figure 4 that the inviscid pres- 
sure distribution is a good approximation to these data, except in 
the vicinity of the shock— wave— boundary- layer interaction, and at the 
trailing edge of the airfoil. 

5.2 Time-Histories of Pressure Coefficients 

Starting from an initial field generated from the inviscid 
small -disturbanGe i solution ..'from' TSFOIL , the zero incidence 
case was run 1100 computational cycles to a characteristic time* ZT 
of 3.62.; The 1100 cycles required 28 minutes on the GDC 76OO computer. 
Pressure— time histories indicated that the iield was near steady- 

state at this characteristic time. 

In order to demonstrate that the airfoil flow field is approaching 
a steady-state, time histories of the pressure coefficient were monitored 
at five points along the surface of the airfoil. The pressure coeffi- 
cient is plotted as a function of the characteristic time T in 
Figure 5. It is seen from Figure 5 that at a characteristic time C 

* Unit characteristic time corresponds to the time it takes a free 
stream particle to travel one chord length. 


of 3.62 the time-histories are nearly horizontal. Although time- 
histories were not generated in other regions of the flow field, it 
is believed that at r = 3.62 the remainder of. the field is also near 
steady-state . ■ 

5.3 Computational Time Reduction Factor 

The STOKES Computer code employed to solve Problem 101.0 embodies 
three computational time reduction methods, namely. 

1. Switched axis 

2. Vector-Do-Loops 

3. Multi -Regional Timesteps 

It is the purpose of this section to determine the speed-up factor 
afforded by each method, and the total speed-up factor. 

let us consider the multi-regional timestep levels (MRT) fiirst. 

Table I presents a tabulation of the timestep levels associated with 
the streamline-like-lines of Figure 3 at a characteristic time C 

of .6600. 

TABLE I 

Tabulation of Timestep Levels 
at a Characteristic Time t of .6606 

Timestep 
us 

.775655 
1.55131 

3.10262 
6.20524 

hence , the computation 
computation. Each streamline- 


K-Line 

Range 


Timestep 

level 


30-33 

25-29 

16-24 

1-15 


1 

2 

3 

4 


There are four timestep levels in Table I; 
reauires eight microsteps per macrostep of 


Number of 
Microsteps 
per Macrostep 

8 

4 

2 

1 


30 


X 


like-line is numbered with the index K, where K =1 corresponds to the 
lower lateral boundary of Figure 3 . It is seen from Table I that four 
K-lines are in Level 1, five K-llnes are in Level 2, nine K-lines are 
in Level 3, and fifteen K-lines are in level 4. On this basis a speed- 
up factor, can be calculated, which is the ratio of the number 

of point-steps that occur if the smallest timestep governs the compu- 
tation to the number of point-steps that actually occur due to MET. 


^MET 


(33)(130)(8) 

130 ( (4)(8) + (5)(4) + ( 9 )( 2 ) + 13J 


= 3A ^ S8 (25) 


The speed-up factor computed above is not the maximum that could be 
achieved for Problem 101.0. The potential-like-lines limit the hxgher 
timestep levels. A revised mesh and a potential-line dropper in the 
STOKES code should give another factor of two. 

The switched axes of STOKES gives a factor of 1.5 (Ref. 14 ) and 
the Vector-Do-Loops employed provides an additional factor of 1.5 
(Eef. l4). Therefore, the total speed-up factor is 


s = ( 1 . 5 ) (1.5) (3.1058) s: 7.00 


( 26 ) 


As was discussed in Section 5.2, the STOKES code was run through 
1100 macrocycles to reach a near steady-state at r=3.82. The 1100 
macrocycles required 28 minutes on the Ames CDC ?600 computer. Based 
on Equation (26), an unswitched constant timestep version of STOKES 
would solve Problem 101.0 in three hours and fifteen minutes on the 
GDC 7600 computer. 

5.4 Velocity Vectors of Steady Flow Field 

Figures 6 and 7 present velocity vector plots of the steady flow 
field about the NACA 64A010 airfoil. The vectors are proportional to 
the particle speeds and emanate from the mesh positions of Figure 3- 
Figure 6 shows the leading edge flow field and the field in the 
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vicinity of the shock wave, while Figure 7 shows the trailing edge 
flow field. 

A deceleration along the dividing streamline of the airfoil is 
shown in the flow field of Figure 6. Furthermore, suhseq^uent expansion 
about the forward portion of the airfoil is also indicated in Figure 6. 
The boundary layer must be thin on the forward portion of the airfoil 
since it is not descemible in the vector plot. 

Aft of the shock wave, which occurs at x =.?4 ft, a thickening of 
the boundary layer is indicated in Figure 6, A boundary layer flow can 
be clearly seen at the axial station x=1.01 ft. 

A recirculation region is clearly indicated in the neighborhood 
of the airfoil trailing edge in Figure ?. This recirculation region 
is only calculable through the Navier-Stokes equations; such a region is 
not within the scope of inviscid methods which employ a Kutta condition. 

The calculated surface pressures are compared to measurements 
in the next section. 

5»5 Surface Pressui^ Distribution Comparisons 

Surface pressure distributions obtained from the STOKES numerical 
solution, the TSFQIl inviscid solution (Ref. 3 )» s^nd the measurements 
of Gross and Steinle (Ref. 34 ) are shown in Figure 8. It is seen that 
the STOKES pressure distribution is nearly Identical to these data, 
except in the vicinity of the shock and at the trailing edge. The 
calculated shock position, i.e., x/c == . 475 , appear correct; however, 
the shock transition is smeared out relative to these data; it is believed 
that the finer mesh computation now in progress will produce a sharper 
shock transition. The TSFOIl prediction produces a shock location 
aft of the experimental value, i.e., x/c = .525, but transition 
is sharper. Finally, at the airfoil trailing edge the TSFOIl prediction 


shows an increase in pressure, while hoth the STOKES and experimental 
results indicate no such increase. The flattening of the trailing edge 
pressure distribution is due to viscous effects. Differences between 
the calculated trailing edge pressures and corresponding data may be 
attributed to an early numerical separation caused by the absence of 
tebulence. It is anticipa-ted that the turbulent computation now in 
progress will greatly narrow differences between the calculated and 
measured pressure fields. 

5.6 Mach Number Contours 

Calculated local Mach number contours are shown in Figure 9. 

The freestream flow at a Mach number of .80 decelerates to a stagna- 
tion point at the leading edge of the system, then accelerates past 
sonic flow to Mach 1.08, shocks down to subsonic conditions, and 
decelerates subsonically throughout the remainder of the airfoil chord. 
Furthermore , a low speed recirculation region exists in the neighbor- 
hood of the trailing edge of the airfoil. 

The gradual transition through the shock wave is indicated in the 
Mach number contours of Figure 9- The airfoil accelerates to Mach 1.08 
at X = .74 ft, then gradually goes through shock transition to 
X = ,94 ft where subsoniq flow occurs. This, smearing of the shock 
will be reduced in the finer mesh calculation now in progress. 


SECTION 6 


RESULTS OF LIFTING CALCULATION 

Computations for the NAGA 64A010. airfoil at two-degrees angle 
of attack were conducted with the turbulence and transition models 
operational. The turbulent lifting airfoil problem was designated 
"Problem 102.0". During the process of running Problem 102.0, it 
was found that an integrated inviscid/viscous calculation provided 
more realistic results than a complete viscous calculation. The 
principal reason for this was that the finite difference mesh employed 
was too coarse to define the thin boundary at the airfoil leading 
edge. The integrated inviscid/viscous calculation was designated as 
"Problem 102. i". 

The lifting airfoil numerical results of Problem 102.0 are presented 
in the following formats 

1. Finite Difference Mesh 

2. Initial Conditions 

3. Computational Time Reduction Factor 

4. General Flow Field Structure 

5. Calculated Lift and Drag 

After the results of Problem 102.0 are described a discussion of 
Problem 102.1 follows. 

6 . 1 Finite Difference Mesh 

The mesh employed to solve Problem 102.0 is presented in Figure 10. 
This mesh is comprised of the intersection of 68 streamline-like-lines 
ajid 130 pbtential-like-lines. The streamline-like-line spacing is 
designed to provide about six points in the leeward boundary layer 
aft of the airfoil shock. 

The mesh of Figure 10 was developed from the half-mesh employed 
for the NACA 64A010 airfoil at zero angle of attack (Figure 3 ) . This 
half-mesh was reflected to produce a symmetric mesh about the full 


64A010 airfoil at zero angle of attack. The airfoil and asso- 
ciated mesh were then rotated through two degrees to produce a finite 
difference mesh at angle of attack. Horizontal lines were then imposed 
at the lateral, "boundaries of the system and the potential-like -lines 
were forced to terminate at these "boundaries . 

In the neigh"borhood of the airfoil leading edge the tur"bulent "boun 
(5.ary layer thickness is much smaller than the thickness of the first 
layer of zones adjacent to the airfoil surface . Therefore ^ the mesh 
of Figure 10 is considered as a medium mesh. It is designed to 
provide a practical running time with sufficient accuracy to define 
the important airfoil fluid mechanical phenomena . 


6.2 Initial Conditions 

The initial field for the solution of Pro"blem 102.0 was 
generated from the inviscld small disturhance solution from TSFOIl 
(Ref. 3) for the 64A010 at M = .8 and 2 degrees angle of attack, 
and interpolated onto the mesh of Figure 10 . 

Initial pressure distri"butions on the leeward and windward sides 
of the airfoil are shown in Figure 11. It is seen from this figure 
that the leeward flow undergoes a very rapid expansion at the airfoil 
leading edge -to supersonic flow, and then is shocked to subsonic 
conditions at a station along the airfoil of approximately of 
the airfoil chord. 

The inviscid lift and drag coefficients are compared to inter- 
polated measurements of Gross and Steinle (Ref. 3^) in Table II. 
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TABLE II 


Comparison of Lift and Drag Coefficients for the NACA 64A010 Airfoil 
Method of Generation Drag Coefficient Lift Coefficient 

2 Ql 

.0413 .8709 

.047 .415 


Murman Inviscid 
Solution (Ref . 3 ) 

Measured Value (Ref, 3^) 


It is seen from Table II that the inviscid theory overestimates the 
lift and underestimates the drag by a large factor. If the Inviscid 
lift prediction were closer to measurements and the flow were steady, 
the Murman far field solution could be imposed on the upstream and 
lateral boundaries of the mesh and the problem run in this way. Due 
to the large discrepancy in lift and the possibility of unsteady flow, 
the boundary conditions described in Section 4.2 were employed on these 
boundaries . 

6.3 Computational Time Reduction Factor 

Starting from the initial field of Section 6.2, the two degree 
incidence case was run 2615 raacrocycles to a characteristic time r 
of 9 •05* The 2615 cycles required 1.88 hours on the GDC 76 OO computer 
at "Ames Research Center. Since Multi-Regional Timestep (MRT) logic 
is employed, each macrocycle is comprised of many microcycles. There- 
fore, the 1.88 hours of computational time is divided according to the 
timestep levels employed in the calculation. The timestep levels 
employed for Problem 102 . 0 are presented in Table III , 
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V 


Macrocvcles 

965 

1650 


TABIE III 

Computer Time Division 
Timesteu Leve'ls 


4 

3 


Computer Time 
( hrs ) 

.857 

1.023 


2615 


1.88 


The results of Table HI indicate that the first 965 macrocycles 
occurred at four times tep levels, i.e., 8 microsteps per macros tep, 
and the remaining I65O macrocycles occurred at three tlmestep levels, 
i.e., 4 microsteps per macros tep. Furthermore, the first 965 cycles 
required ^857 hours of CDC 76OO time, while the remaining I65O cycles 
took 1.023 hours. 

Based on the results of TABIE III the speed-up factor afforded 
by the MRT logic can be evaluated in a manner similar to that of 
Section 5 - 3 * Through the first 965 macrocycles four timestep levels 
were employed. The MRT speed-up factor, Sj^^, that results can be 
determined from Table iv below. 

TABLE IV 

Tabulation of Timestep levels 
for First 965 Macrocycles 


Timestep level 


layere of Zones 
in that level 

12 

15 

28 

12 


1 J 
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130 ( (12)(8) + (28)(2) +12 ) 

The switched axis gives a factor of 1.5 (Ref, ., 13)1 and the Vector-Do- 
Loops provide another factor of 1.5 (Ref. I 3 ) ; therefore, the total 
speed-up factor over the first 9^5 macrocycles is 


= (i.5)(i.5)(2.392) = 5.38 

Through the remaining I 65 O cycles three timestep levels were used. 

'^he MRT speed-up factor , , can be determined from Table V below . 


TABLE V 

Tabulation of Timestep levels for 
Remaining I 65 O Macrosteps 


Timestep Level 

1 

2 

3 


Layers of Zones in the Level 
14 

19 

34 


MRT 


(67)(4)(130 


130 ( (14)(4) + (19)(2) + 34) 


= 2.0938 


Accounting for the switched axis and Vector-Do -loops yields a total 
speed-up factor of 

^650 - 

The overall speed-up, S, is the cycle weighted average of and 


'1650 ’ 


X • 0 • 9 


s = 4.96 

Therefore, in the absence of a switched axis code, Vector-Do-Loops, 
and MRT, Problem 102.0 would require 9.3 hours on the Ames GDC 76 OO 
computer. 
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As was discussed previously, the timestep levels are limited 
hy the potential-lifce-lihe spacing of the finite difference mesh in the 
vicinity of the airfoil leading edge, not the streamline spacing. 
Therefore, zone dropper logic, which drops potential-like-lines of the 
finite difference mesh away from the immediate neighborhood of the 
airfoil leading edge, would increase the number of timestep levels 
from four to five and decrease the number of computations req_uired. 

It is believed an additional speed-up factor of two would result; 
hence, lifting airfoil problems could be solved within an hour. 

6A General Flow Field Structure 

In this section velocity vector plots are shown which illustrate 
the sequence of events as the viscous flow field develops from the 
inviscid initial conditions. Figures 12 through 1? present the velo- 
city field of Problem 102.0 at various characteristic times ranging 
from zero to 8.62. 

Figures 12 to 14 indicate upstream movement of the leeward airfoil 
shock from its initial position to its farthest upstream location. 

Figure 12 shows the initial inviscid field with a leeward shock clearly 
indicated at the 85 percent chord station. At a characteristic time 
r of 2.10 (Figure 13) , the leeward shock has moved to about' the 40 
percent chord station with a separated flow trailing the shock. In 
fact a clockwise vortex is seen about to shed at the trailing edge. 

A large clockwise vortex is shedding from the leeward side in Figure 
(1^3.75). Furthermore, the leeward shock has moved to the 30 per- 
cent chord station, which is approximately its farthest upstream 
location. 

Downstream motion of the shock to a near equilibrium position is 
indicated in the velocity vector plots of Figures 15 through 1?. 

Figure 15 shows the velocity field at a characteristic time r of 5.44. 
At this characteristic time the leeward shock is at the 35 percent 


chord station. The shock has moved to approximately the 50 percent 
chord station at 'C‘=8.27 (Figure 16 );. Figure 1? afC=8.62 also indicates 
a leeward shock position at about the 50 percent chord station. 

The main points of Figures 15 through 1? are that (a) the leeward 
shock seems to arrive at a near-equi librium position, (b) vortex 
shedding appeairs to continue throughout the calculation, and (c) the 
wake of the airfoil assumes a sinusoidal pattern. 

6.5 Calculated lift and Drag 

The main result of the velocity vector plots is that the flow at 
the airfoil trailing edge and in its wake is unsteady. To obtain a 
more quantitative description of this unsteadiness time-histories of 
the lift and drag coefficients are examined. 

The lift coefficient time history is presented in Figure 18. 

Starting from the inviscid lift coefficient of ,8?09 the lift coeffi- 
cient generally decreases in a transient way until a characteristic time 
r* of 7 . 6 . For characteristic’ times greater than 7*6, the lift-time- 
history is periodic , . with a complete period of oscillation shown in 
Figure I 8 . According to inviscid small disturbance theory (Equation 
24-), the circulation is. proportional to lift; hence, the airfoil 
sheds vorticity in a transient way until T of 7-6 a-nd then periodic 
shedding occurs. From Figure 18, the period of oscillation is about 
1 . 4 in characteristic time units , ■ 

A corresponding time history of the drag coefficient is shown in 
Figure I 9 . Starting from the inviscid drag coefficient of .04l3, 
the drag coefficient increases to about .08 at a characteristic time 
C of 3 * 0 » then decreases to about .03 afC= 6 . 0 , and finally becomes 
periodic after a characteristic time of 7*6. The period of oscillation 
of the drag coefficient appears to be similar to that of the lift 
coefficient. 

The lift and drag coefficients of Figure 18 and I 9 , respectively , 
are time-averaged over the period of oscillation. The time-averaged 


lift coefficient is .137 and the time -averaged drag coefficient 
Gj^ is .038. A comparison of these values vfith the corresponding 
interpolated measurements of Table II indicates that the calculated 
values are low. The calculated drag coefficient is in the ballpark, 
but smaller than measured, while the lift coefficient is far less than 
measured . 

The numerical solution was examined to determine the reasons for 
the low lift and drag predictions. The principal reason found was 
that the finite difference mesh in the vicinity of the leading edge 
of the airfoil was too coarse to define the thin boundary layer there. 
The thin leading edge boundary layer is. contained within the first 
layer of zones adjacent to the airfoil on both its leeward and wind- 
ward sides. As a result the high suction pressures near the nose were 
not obtc.ined numerically and the supersonic region was much sitaller than 
that observed experimentally. 

A thin leading edge boundary layer, whose thickness is smaller 
than the width of the first layer of zones of the finite difference 
mesh, has the following two numerical effects; 

1. Calculated airfoil leading edge shear stresses are 
smaller than actual shear stresses; therefore, the 
predicted drag' is less than the actual drag, 

2. Calculated average Mach numbers in the first layer 
of zones about the airfoil leading edge are less 
than they should be. These lower Mach numbers 
permit the periodic flow aft of the leeward shock 
to erode away the expansion region on the leeside . 

An alternate way of explaining the smaller calculated supersonxc 
region is to introduce the concept of displacement thickness^. A 
coarse finite difference mesh at the airfoil leading edge ,xs similar 
to adding a large displacement thickness to the airfoil which inhibits 

the leewax-d expansion. 

* The authors are indebted to Dr, Gary T. Chapman of NASA Ames for 
suggesting this concept. 
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The calculated flow field contains a leeward supersonic region 
and a shock wave; however, the magnitude of the pressure drop due to 
the leeward expansion, and the shock strength are reduced helow the 
actual expansion pressure drop and shock strength. 

There are two methods of correcting this numerical problem. 

1. Be-calculate the flow field with a finer mesh on the 
forward portion of the airfoil solving the full 
Navier-Stokes equations . 

2. Employ another technique, which does not require 

a very fine mesh, to solve for the leeward supersonic 
region, and integrate this solution with the Navier- 
Stokes solution everywhere else. 

It is believed that a fine mesh Navier-Stokes solution is imprac- 
tical at this particular time. The finite difference mesh required is 
about ten times finer than what we are now using. This means that 
the problem will take about 20 hours to solve on the GDC 7^00 computer. 

This computer expenditure makes Navier-Stokes computations about 
lifting transonic airfoils too costly for most applications . 

In this research program a numerical experiment has been conducted 
to investigate the results of integrating another solution with the 
medium mesh Navier-Stokes solution of Problem 102.0. The forthcoming 
section describes this numerical experiment. 

6.6 Integrated Inviscid/Viscous Calculation 

The inviscid steady field from TSFOH , employed 
as the initial conditions, has a much more accurate supersonic region 
than the medium mesh periodic Navier-Stokes solution of Problem 102.0. 

Furthermore, the lesward steady supersonic region computed in Problem 
102.0 during the early stages of motion, i.e., at a characteristic 
time of 1.6, is a more accurate prediction than the TSFOIl field. 
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This is because the airfoil shape is properly accounted for and some 
viscous effects are included. The TSFOIl solution, which employs 
the approximate tangency condition, is interpolated to account for the 
body shape (see Section 5*2). 

The calculated leeward supersonic field of Problem 102,0 at a 
characteristic time of 1.6 was integrated with the Navier-Stokes 
solution everywhere else in two steps. 

1, At characteristic times between ?.6 and 9.05 i*' i'ts 
period of oscillation, the periodic flow field computed 
in Problem 102.0 was patched together with the steady 
leeward supersonic region computed at TT =1.6. 

2. The interface region between these solutions, which 
contained the leeward shock, was recomputed by running 
the ST0X5S code through a characteristic time interval 
of . 025 . 

After the ST0K5S run, an integrated solution resulted which had a 
proper supersonic region, a leeward shock wave-boundary layer inter- 
action, and a separated trailing edge region. The integrated inviscid/ 
viscous calculation was designated as Problem 102.1. 

In the remainder of this section the results of Problem 102.1 are 
presented. The format for presentation is as follows; 

1. General Flow Field Structure 

2. Calculated Lift and Drag 

3. Surface Pressure Distribution Comparisons 

4. Comparison of Mach Number Contours 

5 . Turbulence Characteristic for Airfoil 

6 . 6,1 General Flow Field Structure 

Shedding of vorticity from the lee side of the airfoil is indi- 
cated in the velocity vector plots shown in Figures 20, 21, 22, and 23 . 
Figure 20 presents the velocity field at 'T =7.831, A leeward vortex 
is forming at approximately the 50 percent chord station, while a 


a vortex rotating counter-clockwise is shedding at the airfoil trailing 
edge. The counter-clockwise shed vortex induces clockwise circulation 
about the airfoil? thus, lift is enhanced. At a characteristic time of 
8.204 (Figure 2l), the eye of the leeward vortex of Figure 20 has-moved 
to about the 65 percent chord station and the counter-clockwise shed 
vortex has moved further downstream. Figure 22, at a characteristic 
time of 8.649, shows the leeward vortex of Figure 21 at the airfoil 
trailing edge and in the process of shedding. Furthermore, the sense 
of this vortex is now clockwise. This shed clockwise vortex induces 
a counter-clockwise circulation about the airfoil; hence, the lift 
is reduced. At a characteristic time of 9.052, the velocity field shown 
in Figure 23 no longer has the clockwise leeward vortex in it; however, 
a counter-clockwise vortex appears to be shedding from the trailing 
edge. In fact, the velocity fields of Figures 20 and 23 are somewhat 
similar, suggesting a periodic flow. 

The velocity vector plots of Figures 20 to 23 are similar to the 
Interferogram pictures taken by Finke (Hef. 12) and shown in Figure 1. 
Shedding of vorticity is clearly indicated in both the interferogram 
pictures and velocity vector plots as well as a sinusoidal wake pattern. 

Finally, the alternate shedding of counter-clockwise and clockwise 
ro ta ting vortices in Problem 102.1 indicates that the instantaneous lift 
coefficient should fluctuate with time. Lift and drag are presented in 
the next section. 

6.6.2 Calculated Lift and Drag 

The periodic nature of this flow is quantitatively confirmed in 
calculated time-histories of the lift and drag coefficients shown 
in Figure 24. It is seen from Figure 24 that the instantaneous 

curve has a period of about 1,4 in characteristic time units. 

This period converts to a frequency of 366 hertz, or to a reduced 
frequency* K of 4.49. The drag coefficient time history curve 

* The reduced frequency K is defined as the product of the angular 
frequency and the chord divided by the freestream velocity. 
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ajppeaxB.- to, exhibit the same -frequency. As is shown in the velocity 
vector plots above, during the period of oscillation of the lift 
coefficient alternate shedding of counter-clockwise and clockwise 
vortices occur. 

Independent . experiments on a- symmetric NAC'A 63-j^-012 airfoil at Mach 
number .70, Reynolds number 1.25 x 10^, and two degrees incidence^ 
described in Section 1, show a periodic flow at the trailing edge with 
a reduced frequency K of 2.5 the oscillating shock. The dis- 
crepancy in the reduced frequency between the two cases can be attri- 
buted to differences in airfoil geometry, Reynolds number and Mach number, . 

The calculated instantaneous lift and drag coefficients of Figure 24 
were time averaged and compared to (a) steady lift and drag measurements 
of Gross and Steinle (Ref. 34) and (b) lift and drag predictions of 
small disturbance inviscid theory (Ref. 3 ). The calculated time-averaged 
lift and drag coefficients of Problem 102.1 are compared to correspon- 
ding data in Figure 25. The predicted time-averaged lift coefficient 
is close to the faired experimental curve while the predicted drag 
coefficient is lower. 

A quantitative comparison of predicted lift and drag with inter- 
polated data of Gross and Steinle (Ref. 34 ), and inviscid predictions 

of TSFOIl are shown in Table VI . 


TABLE VI 

Comparison of Lift and Drag. Coefficients for the NAGA 64A010 Airfoil 


Method of Generation 

Drag Coefficient 
°D 

Lift Coefficient 
°L 

Murman Inviscid Solution 
(Small Disturbance) 

.0413 

.8709 

Measured Value 
(Gross arid Steinie, 
Interpolated, TMX -62468) 

.047 

.415 

STOKES Numerical Calculation 
(Problem 102.1) 

.034 

.384 
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It is seen from Table VI that the calculated lift coefficient is within 
seven percent of the mesisured value, while the inviscid prediction is 
greater by a factor of 2.1. Purtheannore, the calculated drag coefficient 
is somewhat lower than measured, while the inviscid prediction is closer 
to the data point. 

The lower calculated drag coefficient is a result of the coarse- 
ness of the finite difference mesh at the leading edge of the airfoil. 

The leeward and windward boundary layers on the NACA 6^lA010 airfoil are 
much thinner than the thickness of the first layer of zones adjacent 
to the airfoil. Thus, the numerically determined shear stress at the 
airfoil leading edge is greatly underestimated, and results in a lower 
drag coefficient prediction. It is believed that for thin boundary 
layers at the airfoil leading edge, boundary layer theory must be inte- 
grated with Navler-Stokes computation everywhere else to produce 
correct drag coefficients. 

6 ,-6. 3 Surface Pressure Distribution Comparisons 

Instantaneous calculated pressure distributions on the leeward 
and windward sides of the NACA 64A010 airfoil are shown in Figures 26 
and 27, respectively. The curves of these- figures correspond to the 
characteristic times 7.831. 8.204, 8.649, and 9.052, which are identi- 
cal to the characteristic times of the velocity vector plots of Figures 
20 to 23. The leeward instantaneous pressure distributions (Figure 26) 
depict the periodic shedding of vortices when compared to the velocity 
vector plots of Figures 20 to 23. 

Let us examine the leeward pressure distributions of Figure 2S. 

Upstream of the 40 percent chord station the pressure field is steady. 

Curve 1, at a characteristic time of 7.831. has a high pressure point 
at the 50 percent chord station followed by a rarefaction whose minimum 
is at the 58 percent chord station. This correlates with the velocity 
vector plot of Figure 20 at r =7.831. Figure 20 shows a stagnation point 
at the 50 percent chord station followed by an expansion about the separated 
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region. At a characteristic time of 8.204 (Curve 2 of Figure 26), 

the high pressure point is at the 60 percent chord station followed 

by a rarefaction; this correlates with the velocity vector plot of 

Figurfi 21, which shows that the stagnation point is at the 60 percent 

chord, station followed by an expansion about the separated region. In 

other words, a separated region has formed at the $0 percent chord 

station atr=?.831 and has moved to the 60 percent chord station atT=6.204. 

Curve 3 of Figure 26 (t= 8.649) shows that the stagnation point has moved 
to the 80 percent chord station with the minimum rarefaction pressure 
at the 90 percent chord station. The velocity vector plot of Figure 
22 (r=8.649) shows that the vortex has moved to the trailing edge of 
the airfoil. Finally, Curve 4 of Figure 26, at a characteristic time 
of 9 . 052 , has no stagnation point followed by a rarefaction. The 
results of Curve 4 indicate that the vortex has shed at r =9. 052. 

The instantaneous windward pressure distributions of Figure 2? 
oscillate near the trailing edge of the airfoil. For percentage chord 
less than 50 percent the flow is steady on the windward side. For 
percentage chord greater than 5° percent oscillations are indicated. 

These windward trailing edge oscillations are a result of the periodic 
vortex shedding process that occurs on the leeward side of the airfoil. 

Instantaneous calculated pressure distributions on the windward 
and leeward sides of the NACA 64A010 airfoil (Figures 26 and 2?) 
were time-averaged over the period of oscillation, A comparison of 
the calculated time-averaged pressures with’ corresponding steady measure- 
ments and with predictions from inviscid small disturbance theory are 
presented in Figure 28. On balance the time-averaged calculated 
pressures compare well with corresponding data, although these data 
are interpolated with angle of attack. The Inviscid theory predicts 
a shock location at the 85 percent chord station while the STOKES 
calculated shock location and the data indicate that the shock is at 
the 45 percent chord station. 
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6,6,4 Comparison of Mach Number Contours 

Mach number contours about the NACA 64A010 airfoil are compared 
in Figures 29 and 30* Figure 29 presents the TSFOIX invlscid 
solution, while Figure 30 represents the STOKES prediction at a charact- 
eristic time of 3=65. The most striking differences between the two 
figures ars the extent of the leeward supersonic region and the locus 
of the leeward shock. The inviscid flow field expands about the leeward 
side of the airfoil to a supersonic Mach number of I.3 and then goes 
through shock transition at the 85 percent chord station. In contrast, 
the viscous field at f*8,65 has a supersonic expansion to Mach 1.2 
and goes through shock transition at about the 45 percent chord station. 

The results of Figures 29 and 30 ars similar with respect to the 
thickness of the shock transition region. Both fields have a leeward 
shock layer approximately 5 percent of a chord thick. The finite 
shock layer of Figure 30 is a result of the STOKES numerical method 
and the finite difference mesh employed in the shock region, if shock 
fitting is not employed, shock transition usually takes place across 
four zones of the finite difference mesh. Hence, the spacing of these 
four zones defines the tnickness of the shock layer. In this particular 
case the layer is approximately 5 percent of a chord thick. A finer 
mesh in the shock region can reduce the thickness of this layer. 

The viscous wake of the airfoil is also indicated in Figure JO, In 
fact, at a station x of about 2.80 the eye of a vortex can be seen 
in the Mach number contour plot. At approximately this station a closed 
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M&ch .80 contour contains a closed Mach .90 contour, indicating a 
decrease in velocity as one moves away from the eye of the vortex. 

Finally, it is pointed out that the Mach number contour plot of 
Figure 30 represents a snapshot of the field at a characteristic time 
of 8 . 65 . The leeward shock actually oscillates and vortex shedding 
occurs. From an examination of the Instantaneous leeward pressure 
distributions of Figure 26, it is seen that the amplitude of the 
oscillation is small,. i.e., about one percent of a chord. This is 
consistent with the data of Finke (Ref. 12 and Section l) for the 
NAGA 63 -012 airfoil at Mach .70 and Reynolds number 1.25 x 10 . 

Plnke found that the shock amplitude was about one percent of a chord 
for 0^=2** and increased to 18^ of a chord at oC=8 . Thus, as the angle 
of attack is reduced towards zero the leeward shock approaches an equi- 
librium position and steady, flow is achieved. 

6 . 6.5 Turbulence Characteristics for Airfoil 

Due to the transient nature of the flow aft of the leeward shock 
wave, it is difficult to assess the validity of the Rose turbulence 
model employed in this numerical investigation. However, numerically 
determined eddy viscosity profiles upstream and downstream of the lee- 
ward shock qualitatively compared with corresponding profiles on 
a flat plate undergoing a shock-wave boundary layer interaction. 

Profiles of the ratio of the eddy viscosity e to the molecular 
viscosity ytt are presented in Figure 31 at two chord stations along 
the leeward side of the NACA 64A010 airfoil. These data correspond 
to a characteristic time ^ of 8 . 65 . The profile on the left hand 
side of the figure, which is at the 31 -^ Percent chord station, is 
upstream of the shock, while the profile on the right hand side of the 
figure, which is at the 6 I.I percent chord station, is downstream of 
the shock. Upstream of^ the shock the profile is nearly linear 
through the laminar sublayer and a portion of the turbulent boundary 


I 




lsiy6f ( constant through the remainder of the boundary layer and then 
rapidly decays as the local frees tream is approached. Downstream of 
the leeward shock the velocity piefile is separated in the vicinity of 
the wall. The separated velocity profile causes the eddy viscosity to 
approach zero near the wall and results in a nearly linear 6 j/A- profile 
in the shear layer above the wall. The linear region is followed by 

a constant value throughout the remainder of the boundary layer and then 
a rapid decay to local freestream conditions. A comparison of the calcu- 
lated profiles upstream and downstream of the shock indicates tnat 

the turbulence intensity increases through the shock transition. 

The (jualitative behavior of the calculated STOKSS profiles is 
similar to corresponding e profiles calculated by Baldwin and Rose 
(Bef. 19 ) on a flat plate. Baldwin and Rose used a relaxation eddy coef- 
ficient model to investigate a shock wave-boundary layer interaction 
on a flat plate at Mach 2.39 and Reynolds number per meter of 5.7 x 10 . 
The Baldwin and Rose profiles upstream of and downstream of the shock 
are presented in the upper right comer of Figure 3 I* -s seen 
from Figure 3 I that the Baldwin and Hose € profiles upstream of and 
downstream of the flat plate shock are similar in shape to corresponding 
NACA 64A010 profiles upstream of and downstream of the leeward 

shock. 
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SEJCTION 7 

CONCLUSIONS AND RECOMMENDATION 

The present investigation has demonstrated that the STOKES computer 
code is an accurate, practical tool for solving viscous, compressible 
lifting and non-lifting airfoil problems at transonic Mach numbers 
suid for Reynolds numbers ranging from the laminar to turbulent regimes . 

Applicability to the turbulent Reynolds number regime was achieved by 
employing the Baldwin-Rose nixing length theory to which relaxation 
along a streamline -like trajectory is incorpo 2 ?ated. 

The numerical method embodied in the STOKES code offers a complete , 
accurate description of the viscous flow field about a transonic air- 
foil, including the important fluid mechanical effects of 

1. shock wave-boundary layer interactions, 

2. transonic buffet and 

3. circulation for the lifting case . 

Furthermore, a Kutta condition, whose formulation for unsteady flew 
is quite uncertaiUjis not reqxdred In this method. 

Accuracy of the numerical method was demonstrated in numerical- 
experimental comparisons for two transonic airfoil cases. The three 
fluid mechanical effects defined above were computed and verified 
experimentally. At zero degrees angle of attack a shock wave-boundary 
layer interaction was successfully computed and verified for the 
NACA 64A010 airfoil. At two degrees Incidence periodic motion (trans- 
onic buffet) was calculated along the aft portion of the airfoil and 
in its wake; this periodic motion has been observed for other airfoil 
geometries at transonic speeds. Although no measurements for the un- 
steadiness of the flow field were made for the NACA 64A010 airfoil at 
these conditions, steady measurements of lift, drag and surface 
static pressures compared favorably with corresponding computed time- 
averaged lift, drag, and surface static pressures. Finally, at two 
degrees incidence the' flow field was computed with far field boundary 
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conditions that simulated the free air case and preserved the airfoil 
circulation throughout the flow field. 

Practicality of the STOKES computer code was demonstrated by the 
Ames GDC 76OO computer time requirements for the two cases solved. 

For the non-lifting case the steady flow field was computed in 28 
minutes, while for the lifting case the periodic flow field was 
computed in 1.88 hours. 

The major criticism of the present method of solution is perhaps 
the accuracy with which the shear stresses are computed at the leading 
edge of the airfoil. Due to the coarseness • of the leading. edge finite 
difference mesh relative to the boundary layer thickness there, the 
calculated leading edge shear stress distribution and hence, the drag, 
was lower than measured. Therefore, it is recommended that in the 
neighborhood of the airfoil leading edge, boundary layer theory be 
integrated with Navier-Stokes computation everywhere else. This will 
produce correct drag coefficients with a minor increase in computational 
time. 
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?lgur« 1. Interferogram pictures of the tine-dependent flow field 
over the NRL 0.1025-0.6750-1.300 airfoil at M -.71, 

amd &C*, “ lxl 06 { pictures depict one full period 
of leeward shock oscillation at Intervals of 700 micro- 
seconds . 
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Figure 6. Velocity vector plot of the flow field on the fonro.irl portion of tfio 
IIACA 6^^A010 airfoil; rj^=.80, x 10^, -0. 
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Ptc«re 7. VelocUy v«lo. pi 7 of the Wllnf e^ro n«. «eW of the f-'.AolO aJrfoll . 
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Figure 9. Local Mach number contoura about the NACA 6UA010 airfoil at zero Incidence t 
AU” * 80 , Re^ ^ X 10^. 
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10. Medium finite difference mesh for NACA 64A010 airfoil at two degrees 
angle of attack. Problem 102.0- mesh comprised of 68 streamlike- 
lines and 130 potential- like-lines . 
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Figure 11. Initial pressure distributions on the leeward and windward surfaces 
of the NACA 64A010 airfoil; M *.8, a- 2°, 4 x 10®. 
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n«u« 15. Velocity yector plot of the^NACA 64A010 etrfoll flo« field at a characteristic tire 
Figure 15. y ^ ^ ^ . 4 , 106 , p„we» 102.0. 













Flgur* 20. Velocity vector plot of the NACA 64A010 airfoil flow field at a characteristic tlwe 
of 7.83t “ .80, -2*' I 6c - - 4 X 106 j Problem 102.1 







Figure 21. Velocity vector plot. of the NACA 64A010 airfoil flow Xleld at a characteristic tli»e 
of 8.20| A^-.80, c<-2'*, ^ X Problem 102.1. 
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Figure 24, Time histories of the lift and drag coefficients for the 

NACA 64A010 airfoil at two degrees angle of attack; Problem 102.1 

integrated inviscid/viscous calculation; M » .80, Re * 4 x 10^. 
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Data of Gross and Stelnle (Interpolated 
TMX62468) 

Time-Averaged Calculated Pressures 
(Problem 102.1) 

Murman Transonic Small Disturbance Theory 
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Figure 28. Comparison of calculated pressures with axperlmental 
data and transonic small disturbance theory; 

M^- .80, Re^ - 4 X 106, o( - 2“. 




Figure 30. Mach number contours about the NACA 64A010 alrfoU at two ^rees angle -of -attack 
and at a characteristic time of 8.651 •80» ^ * ICP i Problem 102.1. 










